An ultracold analogue to star formation: Spontaneous concentration of energy in 

trapped quantum gases 



M. P. StrzyQ and J. R. Anglin 

OPTIMAS Research Center and Fachbereich Physik, 
Technische Universitdt Kaiser slant ern, D-67653 Kaiserslautern, Germany 

Stars form when cold cosmic nebulae spontaneously develop hot spots that steadily intensify until 
they reach fusion temperatures [1 . Without this process, the universe would be dark and dead. 
Yet the spontaneous concentration of heat is exactly what the Second Law of Thermodynamics 
is in most cases supposed to forbid. The formation of protostars has been much discussed [IHS], 
for its consistency with the Second Law depends on a thermodynamical property that is common 
in systems whose strongest force is their own gravity, but otherwise very rare: negative specific 
heat. Negative specific heat turns the world upside down, thermodynamically; it implies that 
entropy increases when energy fiows from lower to higher energy subsystems, opposite to the usual 
direction. Recent experiments have reported negative specific heat in melting atomic clusters [71 18] 
and fragmenting nuclei but these arguably represent transient phenomena outside the proper 
scope of thermodynamics. Here we show that the counter-intuitive thermodynamics of spontaneous 
energy concentration can be studied experimentally with trapped quantum gases, by using optical 
lattice potentials ITO] [TT] to realize weakly coupled arrays of simple dynamical subsystems that 
share the peculiar property of self-gravitating protostars, of having negative micro-canonical specific 
heat. Numerical solution of real-time evolution equations confirms the spontaneous concentration 
of energy in such arrays, with initially dispersed energy condensing quickly into dense 'droplets'. 
We therefore propose laboratory studies of negative specific heat as an elusive but fundamentally 
important aspect of thermodynamics, which may shed fresh light on the general problem of how 
thermodynamics emerges from mechanics. 



The Second Law of Thermodynamics says that entropy 
cannot decrease, which is often summarized by saying 
that heat cannot flow from colder systems to hotter. This 
summary is exact if hotter and colder are interpreted as 
higher and lower temperature, but it is not necessarily 
valid in reference to systems containing more or less en- 
ergy. The rate at which temperature changes with energy 
— the specific heat Cy — is normally positive, so that 
higher temperature implies higher energy, and heat flow 
from high to low temperature redistributes energy more 
evenly throughout aggregate systems, tending towards 
uniform equilibrium. If Cy should be negative, however, 
then higher temperature corresponds to lower energy, 
and increasing entropy requires that systems with less 
energy lose heat to systems with more energy. In aggre- 
gates of many subsystems, negative Cy thus implies in- 
stability toward spontaneous energy concentration. This 
would fit the pattern seen in star formation, but that 
does not prove that the thermodynamic explanation of 
star formation can only be negative Cy- It is debated 
whether or under which circumstances negative Cy is 
really possible (I2ffl6j . 

Some classic thermodynamics texts state flatly that 
systems with negative Cy cannot exist [17], while oth- 
ers discuss the issue at more cautious length [T8 . In 
fact thermodynamics does not prescribe such features of 
any system. It is the discipline of statistical mechanics 
that strives to predict them, by representing the effect 
of complex mechanical interactions in terms of probabil- 



ity distributions, and identifying statistical properties of 
these distributions with thermodynamical quantities like 
Cy. According to the workhorse probability distribution 
of statistical mechanics, the so-called canonical ensemble 
(CE), negative Cy is indeed impossible for any system: 



CE 



> 0, (1) 
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where the averages ( • )ce are taken in the canonical en- 
semble (see Supplementary Material). This theorem does 
not hold, however, in the alternative probability distri- 
bution of the micro-canonical ensemble (/iC), in which 
the system's energy is flxed, rather than fluctuating ran- 
domly as in the CE. It has been argued that the canonical 
ensemble is invalid whenever the micro-canonical ensem- 
ble yields negative Cy [19 . These discrepancies between 
ensembles point to the fact that statistical mechanics re- 
mains a comparatively weak link in the chain of physics. 

It is nonetheless a tremendously important link. Sta- 
tistical mechanics gives physics much of its power, by 
connecting simple models to complex reality, and its pre- 
dictions are amply verifled in many cases. They are based 
on falsiflable assumptions, however, and these can be 
tested today as never before, with experiments on tightly 
controlled mesoscopic systems whose dynamics can be 
followed closely enough for direct comparison with the 
statistical theory. Trapped quantum gases offer an es- 
pecially promising laboratory for such tests. The ob- 
servation of Bose-Einstein condensation was a scientiflc 
landmark because it conflrmed one dramatic prediction 
of statistical mechanics, the cessation of thermal motion 
by a large fraction of gas atoms at a temperature above 
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FIG. 1. (Color online) Two-mode Bose-Hubbard systems 
weakly coupled in a two-dimensional sheet. Atoms oscillate 
between modes within each subsystem at the tunnelling Rabi 
rate Q, as well as between the subsystems at cj <^ O. 



absolute zero. Can trapped quantum gases test another 
dramatic statistical mechanical prediction, spontaneous 
energy concentration as in star formation, by realizing an 
aggregate of systems with negative specific heat? 

They can. The two-mode Bose-Hubbard (BH) model 
with repulsive interactions has been highly studied as a 
model system, representing bosonic particles that quan- 
tum mechanically tunnel back and forth between two po- 
tential wells [20l |2T] . It has also been realized as such in 
experiments [22H25] . In the formalism of so-called second 
quantization, where the canonical operators 2 remove 
atoms from wells 1 and 2, and their conjugate operators 
a\ 2 replace the atoms correspondingly, the Hamiltonian 
of a single two-mode BH system reads 



2 

cr=l 



(2) 



For repulsively interacting trapped bosons, the constants 
> and £ > are determined by the details of the 
trapping potential and of interparticle scattering, respec- 
tively, is chosen equal to the expectation value of 
a\a-^ + a2a2, which represents the total number of par- 
ticles in both wells together, and is conserved by time 
evolution under J^2- 

A system with only two degrees of freedom like this 
is not normally treated statistically, but if we consider 
a large collection of many such systems, weakly coupled 
together as sketched in Fig. [l] then applying statistical 
mechanics to the individual two-mode systems is exactly 
like the standard textbook case of deducing the thermo- 
dynamic properties of a dilute gas from the statistical 
mechanics of a single non-interacting particle. So with 
precisely the same logic, we obtain the quantum statis- 
tical mechanical entropy of two-mode BH in the micro- 
canonical ensemble as 



where Z/\e{E) is the number of eigenstates of H2 with 
eigenvalues in the range {E^E + AE)^ and ks is Boltz- 
mann's constant. Here the ground state energy and en- 
tropy have both been set to zero by shifting all energies 
and entropies uniformly, as one is free to do in thermo- 
dynamics. In the thermodynamic limit of a dense energy 
spectrum, where this definition is normally applied, AE 
must be chosen great enough for Z/\e{E) to be large and 
relatively smooth as a function of E^ but small enough 
for Z/\e{E) to be linearly proportional to AE for all E, 
so that Sae{E) S{E) becomes independent of AE. 

The two-mode BH spectrum becomes dense in the limit 
of large total boson number which is also the regime 
most easily attained in experiments, and in this semiclas- 
sical limit the eigenspectrum is given accurately, except 
for classical orbits too near an unstable fixed point, by 
Bohr-Sommerfeld quantization. This involves construct- 
ing canonical action-angle variables J, (j) for the quan- 
tum system's classical analogue, expressing the classical 
Hamiltonian H2 as E{J)^ and then determining the quan- 
tum eigenenergies as En = E{27rHn) for n a whole num- 
ber. The limit of a dense spectrum is obtained for large 
N, and implies Zae{E) = {dJ/dE)AE/{27TH). For the 
classical analogue of H2 the derivative dJ/d£^ may be 
calculated analytically in closed form (see the Supple- 
mentary Material), yielding the entropy 



S{E) = kB In 



-2 K[k{E/E^,^,e)] 
^eR{E/E^,^,e) 



(4) 



where K{k) is the complete ehiptic integral of the first 
kind, and we define the functions 



(5) 



Sae{E) = kBHZAE{E)/ZAEm 



(3) 



H£, e) = - ^[e-e{R{£,e)-l/ef]/R{£,e) . (6) 

Here £^max(^, ^) denotes the highest energy of the clas- 
sical system for fixed so that £ = E/E^g,^ is the 
normalized energy. The entropy turns out only to de- 
pend on N and Q by depending on £. (In this sense 
the two-mode BH entropy is non-extensive with parti- 
cle number but this should not be surprising, since 
what is expected is extensivity in the number of two- 
mode systems that are coupled together in a large ag- 
gregate system.) The specific heat may then be calcu- 
lated microcanonically in terms of derivatives of S{E): 
Cf = -{dS/dE)y{d^S/dE^). 

Both S and Cf.^ are plotted versus £ for different val- 
ues of s in Fig. pi Note that Cy^ is always negative. For 
6 = 1, Cy^ diverges at £^ = 1, but at this point Bohr- 
Sommerfeld quantization breaks down for any N because 
the maximum energy classical orbit is a dynamically un- 
stable fixed point. For e < 1 and large the statistical 
mechanical result is definite: an array of weakly coupled 
two- mode BH systems, as in Fig. [l] is an aggregate of 
interacting subsystems that all have negative Cy- 
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FIG. 2. Left: Microcanonical Entropy S{E) of the two-mode 
BH system in units of /cb as a function of normalized energy 
for s = 0.2,0.4,0.6 and 0.8 (bottom to top), where ground 
state energy and entropy are both set to zero. Right: Micro- 
canonical heat capacity Cv{E) of a single two- mode system 
as a function of energy for the same values of e (top to bottom 
at the left of the graph) 



Realizing a large array of two- mode BH subsystems, 
each weakly coupled to its neighbours, is certainly pos- 
sible in current laboratories that trap cold bosons in 
so-called optical lattice potentials, induced by standing 
waves of laser light [22] . The weak coupling, which allows 
the subsystems to exchange both atoms and energy with 
their neighbours, is provided by quantum tunnelling, just 
as between the two wells within each two-mode subsys- 
tem, except that a slower frequency <C 1^ is imposed, 
by making the potential barriers higher. 

It is also possible to excite some of these two-mode 
subsystems locally. If the assumptions of statistical me- 
chanics are correct, then thermodynamics implies that 
such initially distributed energy will concentrate, as it 
does in star-producing nebulae, rather than dispersing 
as it does in almost all other cases, because the weakly 
coupled subsystems all have negative specific heat. This 
is a strict test of statistical mechanics, but it is a fair 
test. Predicting the behaviour of complex aggregate sys- 
tems like the array shown in Fig. [l] from thermodynamic 
properties of the isolated components of the aggregate as 
plotted in Fig. |2j is precisely what statistical mechanics 
is supposed to do. Does it pass the test? 

In a basic regime we can confirm that it does. The 
Hamiltonian for the large aggregate system reads 

H = T.^v+T.p (7) 

2 

a—1 

^iJ = "T ^ (4,y<i+l,i + + l + H- C.) . 

a 

in which we recognize as a two-dimensional array of 

H2 forms, with the nonlinearity constant written with U 
rather than £HQ/{2N) because N Nij can now vary 



in time, as the T-j term lets particles tunnel between 
neighbouring subsystems. In the experimentally attain- 
able limit where all the Nij remain large, the complex 
quantum many-body dynamics of H can be well approx- 
imated with classical mean- field theory [26] . Even the 
classical dynamics is chaotic, but it can be integrated 
numerically. The results for a 512-by-512 array of two- 
mode BH systems are shown in Fig. [3j as three 'stills' 
from a video available in the Supplementary Material. 

The non-equilibrium initial state, which includes only 
weak energetic inhomogeneity, relaxes by excitation of 
sound waves, i.e. waves in the number distribution Nij, 
seen in the left column of Fig. [Sj Some of the initial 
energy thus disperses into these complex low-amplitude 
waves; but the majority of the two-mode energy concen- 
trates into bright droplets of maximum excitation, seen in 
the right column of the Figure. The droplet boundaries 
are also visible in the left column as domain- wall-like de- 
pressions in the local particle number, which extend over 
several lattice sites. Their thickness is on the order of 
the characteristic 'healing' length of the mean-field the- 
ory [26 . The regions of energy concentration behave very 
much like droplets, with positive surf ace te nsion. They 
can also move (see the Supplementary video). 

The results shown in Fig. |3] and the Supplementary 
[video] are robust and typical over a wide range of parame- 
ters and initial conditions. The numerical analysis there- 
fore confirms the dramatic statistical mechanical predic- 
tion of spontaneous energy concentration in the semiclas- 
sical limit. This is our first main result. Moreover, we 
are able to follow the formation of 'droplets of heat' as a 
dynamical process. At least in the regime we have anal- 
ysed here, we can thereby 'reverse engineer' the thermo- 
dynamics of energy concentration, in the sense that we 
can identify the dynamical mechanisms underlying the 
phenomenon. 

The weak coupling limit u <^ ft implies a time scale 
separation, such that when two neighbouring subsystems 
exchange energy, the sum of their two action variables J 
is a so-called adiabatic invariant. In a Bohr-Sommerfeld 
semiclassical sense, this implies that quanta of local two- 
mode BH energy are approximately conserved. A re- 
summed Bogoliubov transformation [27l [28] shows the 
same result perturbatively, rather than semi-classically. 
One can then derive a low-frequency effective theory for 
the evolution of atoms and plasmon-like two-mode ex- 
citations, as two species of separately conserved quasi- 
particles. This comes near to reviving the 18th century 
caloric theory of heat [29 , in a quantum mechanical re- 
interpretation. 

Expanding the full Hamiltonian H in terms of these 
new quasi-particles, we find that the leading interactions 
among the plasmon-like energy quasi-particles are at- 
tractive whenever those among the atoms are repulsive 
[271 128] . Within this purely mechanical effective theory, 
the 'heat droplets' predicted by statistical mechanics thus 
appear as structures like the spin domains seen in spinor 
condensates ^30^ , which are understood in terms of purely 
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FIG. 3. (Colour online) Left column: Total particle number in 
each two-mode system Nij at three points in time (t = Or, 92t 
and 697t, where r = 2n/uj). Right column: Normalized en- 
ergy of the two-mode system E / Ema^ at each site at the same 
points in time. System parameters: = /icj^O.l, [7 = 2 
and N — 1024, lattice size: 512x512. We start with uniformly 
distributed atoms and diffuse clouds of excitation defined by 
simulated phase imprinting. These clouds concentrate into 
maximally excited droplets. 



mechanical instabilities. In repulsive BH systems, energy 
concentrates spontaneously because local energy excita- 
tions attract each other. 

This promising insight into the microphysical basis 
of an important thermodynamic phenomenon has only 
been achieved in a simple limit, however. Beyond the 
mean- field or perturbative regimes, the quantum many- 
body theory of large BH arrays with large particle num- 
ber and high excitation is extremely challenging, even 
with the best available computations. Spontaneous en- 
ergy concentration is a robust enough mean- field effect. 



however, that it must surely extend in some form into 
more strongly quantum mechanical regimes. Our sec- 
ond main contribution here is to have shown how much 
can be learned about the mesoscopic interplay between 
quantum mechanics and thermodynamics, by exploring 
the dramatic phenomenon of spontaneous energy concen- 
tration, with quantum gas experiments on coupled Bose- 
Hubbard systems. We can bring the heat that kindles 
stars into the ultracold laboratory. 



SUPPLEMENTARY MATERIAL 

Time evolution video 

The Supplementary video of the evolution of the local 
particle number and the local energy ofa512x512 lattice 
of two-mode BH systems can be found online at 



http : //www . physik . uni-kl . de/ uploads/media/ 
Strzys-Droplet-movie .mp4[ 

For more information cf. the caption of Fig. 3. 



A. Specific heat in the CE 

The expectation value of the energy in thermal equi- 
librium according to the CE is given by 



(^)cE = ^^zexp(-/3^,) 



^exp(-/3^,) 



(8) 



where (3 = l/{k-BT). The specific heat may then imme- 
diately be calculated according to 



(9) 



= k^P\{E-{E)f)c^ 
which is equivalent to equation ([T]). 

B. Entropy and negative specific heat 



It is easy to show that two systems with negative heat 
capacity cannot be in ordinary thermal equilibrium. To 
see this suppose that total entropy S of two subsystems 
with total energy E is given by 



S = Si{Ei) ^ S2{E - El) 



(10) 



In equilibrium we should have maximum entropy for 
given total energy, and we require dS/dEi = to have 
an entropy extremum, and hence the temperatures of the 
subsystems Ti = {dS/dEi)~^ must be equal, Ti = T2. 
But this is only equilibrium if this extremum of the en- 
tropy is a maximum. Therefore we must consider also 
the curvature of the entropy, given by 



dEl 



1 



1 



(11) 
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with Xjd = -T^{d^S/dEf). If both Ci < and C2 < 0, 
we always have d'^S/dEf > 0, and 5* has a minimum at 
Ti = T2. Any entropy-increasing but energy-conserving 
fluctuations will spontaneously bring the system away 
from the equal temperature state, by heating up one of 
the subsystems at the price of cooling the other one down. 
If the double system does settle down to a maximum en- 
tropy state, the two subsystems will be strongly corre- 
lated, with one having most of the energy, and the other 
very little. Such states cannot be described by the CE, 
because its defining feature is that it makes the subsys- 
tems uncorrelated. In an aggregate of many subsystems, 
however, a set of subsystems which happen all to have 
nearly equal energies may admit a CE description locally, 
allowing a canonical definition of temperature that would 
correspond to the local measurements furnished by a suf- 
ficiently small thermometer. 



C. Entropy of a two- mode BH system 

To compute the micro-canonical entropy of two-mode 
BH in the Bohr-Sommerfeld limit, we must construct the 
classical action-angle variables for this system. It is more 
convenient to do this by abandoning the explicit creation 
and destruction operators in favour of angular momen- 
tum according to the Schwinger representation: 













— a^a 




fl /,t~ 
= 2 - 





(12) 



With these definitions (j2| can be recast into 



(13) 



where the last term, depending only on TV, commutes 
with the rest of the Hamiltonian. Since it therefore only 
shifts all energies by a constant, it will be omitted with- 
out loss of generality in the following. 

The classical energy of the two-mode BH system thus 
is equal to 



E 



(14) 



We then introduce new variables R > and 7 through 
the ansatz 



= ^ ( itcos(7) - - 



Nh 



i?sin(7). 




(15) 



(16) 



The ± in the definition of means that we use two 
patches of the i?, 7 co-ordinates to cover the Lx^y^z 
sphere. Note that the requirement that all three of L^^y^z 
remain real constrains the ranges of R and 7. In partic- 
ular, for £ < 1, must lie between — 1 and + 1. 
In terms of the new variables the energy reads 



E 



NHQ- 
4 



1 



E^,^S{R,e)^Eo, (17) 



independent of 7. This shows that is a constant of the 
motion, and so i?, 7 are already a step towards the action- 
angle variables, and only need to be rescaled. In ( p!7| ), 
^max + Eq and £^0 are by definition the maximum and 
minimum values of E{R,£)^ so that E{R^e) G [0,1] by 
construction. Since the maximum and minimum of E are 
attained, respectively, for the minimum and maximum 
values of R^ namely ^1, we we see that the maximum 
excitation energy is £^max = NhQ^ and £^0 = —NhVt/2. 
Solving (17) then for R{£,e), we find 



R{S,e) 



4eS. 



(18) 



To determine the proper rescaling that makes true 
action- angle variables, we now use the canonical equa- 
tions of motion for the R^ 7 variables, which read 



R = 0, j = €Q\l 



1 



R^^-Rcos{j). (19) 



Defining the elliptic parameter 
1 



Vl-{R-e-^y 



(20) 



the equation of motion of the angular variable 7 can be 
integrated to yield an analytic expression for the period 



ri£,e) = 



8k{£,e)K{k{£,e)) _ 4K{k) 



(21) 



of a classical orbit; here K{k) is the complete elliptic 
integral of the first kind. On the other hand, in terms 
of the properly rescaled action-angle variables J, of the 
classical system, the canonical equations must read 



J = 0, 



dE 



27T 



dJ T{R,e) ' 

This implies that for our system we must have 

a/ _ 7^ _ 2K{k) 
dE~2^~ 7ry^eR{£,e)n ' 



(22) 



(23) 



Using this result we may therefore analytically cal- 
culate the number states within AE of E in Bohr- 
Sommerfeld quantization to be 



Z^e{E)-— — 



K{k) AE 



(24) 
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Using the identity i^(0) = 7r/2, we compute the micro- 
canonical entropy with both entropy and energy of the 



ground state set to zero to be 

S{E) = kB In [ZAE{E)/ZAEm 
= kB In 



2K{k) VlTe 



leR 



(25) 



as quoted in our main text, for k{E^e) given by (20). 
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